#include "meshTypeFunctions.hpp"
#include "meshImportExport.hpp"
#include "uniMesh.hpp"
#include "rectMesh.hpp"
#include "structMesh.hpp"
#include "mesh.hpp"
#include "streamAccess.hpp"
#include "topology.hpp"
#include "util.hpp"
#include <vector>
#include <utility>
#include <list>

namespace hms
{

template<>
Mesh uniMeshAs(
	Array2i meshSize,
	Array2s diff,
	Array2s orig,
	bool addBoundaryPatches
){
	Array2i nVerts { meshSize+1 };
	
	Matrix2Xs vertices (2, nVerts.prod() );
	
	for (Index j {0}; j<nVerts.y(); ++j){
		for (Index i {0}; i<nVerts.x(); ++i)
			vertices.col(j*nVerts.x() + i) = orig + ( diff * Array2s{i,j} );
	}

	Index nCells { meshSize.prod() };
	ArrayX<char> nDomainCellVerts { nCells };
	Array4Xi domainCellVerts { 4, nCells };

	nDomainCellVerts = 4;

	Index nx { meshSize.x() }, ny { meshSize.y() };

	for ( Index j {0}; j < ny; ++j ){
		for ( Index i {0}; i < nx; ++i ){
			/* put vertex indices counterclockwise */
			domainCellVerts.col( j*nx + i ) = Array4i{
				j*(nx+1) + i,         // bottom left
				j*(nx+1) + i + 1,     // bottom right
				(j+1)*(nx+1) + i + 1, // top right
				(j+1)*(nx+1) + i      // top left
			};
		}
	}

	Mesh mesh {
		std::move(vertices),
		std::move(nDomainCellVerts),
		std::move(domainCellVerts)
	};

	if (addBoundaryPatches){
		Index tl { (nx+1)*ny }, tr { (nx+1)*(ny+1) - 1 };
		Array3i
			bot {0, 1, nx},
			lft {0, nx+1, tl},
			rgt {nx, 2*nx + 1, tr},
			top {tl, tl+1, tr};
		
		Array3Xi edgesAtBoundaryVerts { edgesAtBoundaryVertices(mesh) };

		auto addPatch = [&](std::string&& name, const Array3i& verts ){
			mesh.addBoundaryPatch(
				name,
				verticesToBoundarySegment( mesh, verts, edgesAtBoundaryVerts )
			);
		};

		addPatch("left", lft);
		addPatch("right", rgt);
		addPatch("top", top);
		addPatch("bottom", bot);
	}

	return mesh;
}

template<>
UniMesh uniMeshAs<UniMesh>(
	Array2i nCells,
	Array2s diff,
	Array2s orig,
	bool addBoundaryPatches
){
	UniMesh mesh { nCells, diff, orig };
	if (addBoundaryPatches)
		addStandardBoundaryPatches(mesh);
	return mesh;
}

template<>
RectMesh uniMeshAs<RectMesh>(
	Array2i nCells,
	Array2s diff,
	Array2s orig,
	bool addBoundaryPatches
){
	VectorXs xCoords (nCells.x()+1);
	VectorXs yCoords (nCells.y()+1);
	
	for (Index i{0}; i<xCoords.size(); ++i){
		xCoords(i) = orig.x() + i*diff.x();
	}
	for (Index i{0}; i<yCoords.size(); ++i){
		yCoords(i) = orig.y() + i*diff.y();
	}
	RectMesh mesh {std::move(xCoords), std::move(yCoords), false};

	if (addBoundaryPatches)
		addStandardBoundaryPatches(mesh);
	return mesh;
}

template<>
StructMesh uniMeshAs<StructMesh>(
	Array2i nCells,
	Array2s diff,
	Array2s orig,
	bool addBoundaryPatches
){
	Matrix2Xs vertices (2, (nCells+1).prod() );

	for (Index j{0}; j<nCells.y()+1; ++j){
		for (Index i{0}; i<nCells.x()+1; ++i){
			vertices.col(j*(nCells.x()+1)+i) =
				Array2s{i,j}*diff + orig;
		}
	}
	StructMesh mesh { std::move(vertices), std::move(nCells) };

	if (addBoundaryPatches)
		addStandardBoundaryPatches(mesh);
	return mesh;
}




auto uniqueCoordinates( const Ref<const Matrix2Xs>& vertices) ->
	std::pair<VectorXs, VectorXs>
{
	std::vector<scalar> xCoords;
	std::vector<scalar> yCoords;

	auto searchAndInsert = [](std::vector<scalar>& v, scalar s) -> void {
		auto it {std::upper_bound(v.begin(),v.end(), s)};
		/* if no entry greater or equal was found */
		if ( it == v.end() ){
			if ( v.empty() || !approxEqual( *(it-1), s ))
				v.push_back( s );
			return;
		}
		if (
			   approxEqual(  *it   , s )
			|| approxEqual( *(it-1), s )
		)
			return;
		else
			v.insert( it, s );
	};

	for ( const Ref<const Matrix2Xs>::ConstColXpr vertex : vertices.colwise() ){
		searchAndInsert(xCoords, vertex[xdim]);
		searchAndInsert(yCoords, vertex[ydim]);
	}

	return {
		Eigen::Map<VectorXs>(&xCoords[0], xCoords.size() ),
		Eigen::Map<VectorXs>(&yCoords[0], yCoords.size() )
	};
}


bool isQuadrilateral( const Mesh& mesh ){
	return (mesh.cells().nDomainCellEdges().array() == 4).all();
}

/**
 * @brief Determines whether a mesh is structured, i.e. has regular neighbouring
 * relations and could be represented by a StructMesh object
 * 
 * @param mesh
 * @param edgesAtVertex the return of edgesAtVertices, if available
 * @return std::optional<IdxVector> If the passed in Mesh is structured, the
 * returned std::vector<Index> holds the vertex indices of the four corners of
 * the Mesh. If it is not structured, the return object does not hold a value
 * and thus converts to boolean false.
 * 
 * Uses the number of edges connected to each vertex to make the determination.
 * Each inner vertex must have four connecting edges, each boundary vertex must
 * have three, each corner two.
 */
auto isStructured(
	const Mesh& mesh,
	const std::vector<IdxVector>& edgesAtVertex
) ->
	std::optional<IdxVector>
{
	if ( !isQuadrilateral(mesh) ) return std::nullopt;
	
	Index vertex {0};
	IdxVector corners;

	for ( const auto& edges :
		( edgesAtVertex.empty() ? edgesAtVertices(mesh) : edgesAtVertex)
	){
		bool innerVertex { true };
		
		for ( Index edge : edges ){
			if ( mesh.edge(edge).isBoundary() ){
				innerVertex = false;
				break;
			}
		}
		
		if (innerVertex) {
			if ( edges.size() != 4)
				return std::nullopt;
		}
		else if (edges.size() > 3)
			return std::nullopt;
		else if (edges.size() == 2) {
			corners.push_back(vertex);
		}
		++vertex;
	}
	if (corners.size() != 4)
		return std::nullopt;
	else
		return corners;
}

/* this overload requires the unique x- and y-coordinates to already be
 * determined. This is the one called by optimiseMesh() */
bool isRectangular(
	Index nDomainCells,
	const Ref<const VectorXs>& xCoords,
	const Ref<const VectorXs>& yCoords
){
	return (xCoords.size()-1)*(yCoords.size()-1) == nDomainCells;
}

bool isUniform(
	const Ref<const VectorXs>& xCoords,
	const Ref<const VectorXs>& yCoords
){
	assert( xCoords.size() > 1 && yCoords.size() > 1 );

	scalar dx { xCoords[1] - xCoords[0] };
	scalar dy { yCoords[1] - yCoords[0] };

	auto uniformLength = [](const Ref<const VectorXs>& v, scalar s) -> bool {
		for ( Index i{1}; i<v.size()-1; ++i){
			if ( !approxEqual( v[i+1] - v[i], s ) )
				return false;
		}
		return true;
	};

	if ( uniformLength(xCoords, dx) && uniformLength(yCoords, dy) )
		return true;
	else
		return false;
}



StructMesh meshToStructMesh(
	Mesh&& mesh,
	const std::vector<IdxVector>& edgesAtVertex,
	const IdxVector& corners
){
	/* Two tasks must be accomplished to convert a Mesh to a StructMesh:
	 * 1. reorder the vertices to conform to the standard which StructMesh
	 *    relies on (left->right, bottom->top)
	 * 2. add equivalent boundary patches
	 * */

	/* a StructMesh is defined by its vertex coordinates. The number is the same
	 * as in the unstructured mesh, but the order can be different */
	/* This vector is for the index translation of vertices from Mesh to 
	 * StructMesh, e.g. at position 0, the first vertex of the to-be-created 
	 * StructMesh, it holds the index of the corresponding vertex within the
	 * Mesh */
	ArrayXi vertices ( mesh.nVertices() );

	/* the starting point is the bottom left corner */
	vertices(0) = findBottomLeft( mesh.vertices()(hms::all, corners) );

	/* we want to go right, so the next vertex is the one further to the right*/
	Vector2i attached {
		otherVertex(mesh, edgesAtVertex[vertices(0)][0], vertices(0)),
		otherVertex(mesh, edgesAtVertex[vertices(0)][1], vertices(0))
	};
	vertices(1) = (
		mesh.vertex( attached[0] )[xdim] >
		mesh.vertex( attached[1] )[xdim] ) ?
		attached[0] :
		attached[1];

	Index nCellsX {1};
	/* the number of cells in x- and y-direction is as of yet unknown */
	while (true){
		if ( edgesAtVertex[vertices(nCellsX)].size() == 2 ) break;
		for ( Index edge : edgesAtVertex[vertices(nCellsX)] ){
			/* we want to go along the boundary to the next corner, so the first
			 * condition is not to go backwards, and the second condition is to
			 * stay along the boundary */
			Index candidate { otherVertex(mesh, edge, vertices(nCellsX)) };
			if (
				candidate != vertices(nCellsX-1) &&
				edgesAtVertex[candidate].size() < 4
			){
				/* not a mistake. write to next position */
				vertices(++nCellsX) = candidate;
				break;
			}
		}
	}
	Index nCellsY { mesh.nDomainCells() / nCellsX };
	assert( mesh.nDomainCells() % nCellsX == 0 );

	/* parallelisable only within one row. Chunk size as nCellsX-1 / nproc,
	 * then barrier */
	for (Index i = 0; i<mesh.nVertices() - (nCellsX+1); ++i){
		for ( Index edge : edgesAtVertex[vertices(i)] ){
			/* In the previous loop, the first row of vertices in topological
			 * x-direction was established. To get to the next row, we must go
			 * in topological y-direction ("up").
			 * Candidates are all vertices connected by an edge to the current
			 * row. The correct one is the one in the topological upwards-
			 * direction. The if-clause is to exclude all other candidates, i.e.
			 * those which are in the current row and the one below.
			 * While in the first row, no row below exists and the very first 
			 * vertex does not have a left neighbour, so going out of bounds
			 * in these cases must be guarded against. */
			Index candidate { otherVertex(mesh, edge, vertices(i)) };
			if (
				candidate != vertices(i+1) &&
				(i >= 1         ? (candidate != vertices(i-1)) : true ) &&
				(i >= nCellsX+1 ? (candidate != vertices(i-nCellsX-1)) : true )
			){
				vertices(i+nCellsX+1) = candidate;
				break;
			}
		}
	}
	/* the vertex matrix for the new mesh (via indexed view) */
	Matrix2Xs structVertices (2, mesh.nVertices() );
	structVertices = mesh.vertices()(hms::all, vertices);

	StructMesh structMesh { std::move(structVertices), nCellsX, nCellsY };

	/* for boundary patch translation, the opposite vertex translation direction
	 * is required. verticesBackwards holds at position i the vertex index of
	 * the mesh's ith vertex within the structMesh, e.g. if the vertex with the 
	 * coordinates (1,2) is at index 3 within the Mesh and at index 4 within
	 * the StructMesh, then verticesBackwards[3] = 4 */
	ArrayXi verticesBackwards ( vertices.size() );
	for (Index i = 0; i<vertices.size(); ++i){
		verticesBackwards(vertices(i)) = i;
	}

	for ( const BoundaryPatch& patch : mesh.boundaryPatches() ){
		Array3i vertsMesh { boundarySegmentToVertices( mesh, patch.edges() ) };
		Array3i vertsStructMesh { verticesBackwards( vertsMesh ) };
		addBoundaryPatch( structMesh, patch.name(), vertsStructMesh );
	}
	return structMesh;
}



/**
 * @brief Checks whether a mesh can be represented by a lower-complexity mesh
 * type and, if so, performs the conversion
 * 
 * @param meshVar a std::variant of all mesh types, with the active one being
 * the current representation.
 */
bool optimiseMeshType( MeshVariant& meshVar ){

	bool optimisationOccurred {false};
	/* lambdas for reusage on different mesh types */
	/* furthest optimisation, rectangular->uniform */
	auto rectOpt = [&]( auto& rectMesh, auto& xCoords, auto& yCoords ) -> bool {
		if ( isUniform(xCoords, yCoords) ){
			UniMesh uniMesh {
				xCoords.size()-1,
				yCoords.size()-1,
				xCoords[1] - xCoords[0],
				yCoords[1] - yCoords[0],
				Vector2s{xCoords[0], yCoords[0]}
			};
			rectMesh.moveBoundaryPatchesTo(uniMesh);
			meshVar = std::move(uniMesh);

			optimisationOccurred = true;
			return true;
		}
		return false;
	};
	/* structured -> rectangular (-> uniform) */
	auto structOpt = [&]( StructMesh& structMesh ) -> void {
		auto [xCoords,yCoords] { uniqueCoordinates( structMesh.vertices() ) };
		if ( isRectangular(
			structMesh.nDomainCells(),
			xCoords,
			yCoords
		)){
			/* check whether mesh is uniform, else convert to rectangular */
			if ( !rectOpt( structMesh, xCoords, yCoords ) ){
				RectMesh rectMesh {
					std::move(xCoords), std::move(yCoords)
				};
				structMesh.moveBoundaryPatchesTo(rectMesh);
				meshVar = std::move(rectMesh);

				optimisationOccurred = true;
				return;
			}
		}
	};

	/* here the actual checks and optimisation calls begin */
	if ( std::holds_alternative<Mesh>(meshVar) ){
		Mesh& mesh { std::get<Mesh>(meshVar) };
		std::vector<IdxVector> edgesAtVertex { edgesAtVertices(mesh) };
		/* perform checks from least to most optimised mesh type, as they are
		 * only valid in that order */
		if ( auto corners = isStructured(mesh, edgesAtVertex) ){
			/* convert to StructMesh first, because the indexing is required
			 * to maintain boundary patches */
			StructMesh& structMesh = meshVar.emplace<StructMesh>(
				meshToStructMesh( std::move(mesh), edgesAtVertex, *corners )
			);
			structOpt( structMesh );
			optimisationOccurred = true;
		}
	}
	else if ( std::holds_alternative<StructMesh>(meshVar) ) {
		structOpt( std::get<StructMesh>(meshVar) );
	}
	else if ( std::holds_alternative<RectMesh>(meshVar) ) {
		RectMesh& mesh { std::get<RectMesh>(meshVar) };
		/* return value can be ignored */
		rectOpt( mesh, mesh.xCoords(), mesh.yCoords() );
	}
	/* if it's of type UniMesh, it already has the most memory-efficient type */
	return optimisationOccurred;
}

} // namespace hms
